Real space tests of the statistical isotropy and Gaussianity of the WMAP CMB data 
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We introduce and analyze a method for testing statistical isotropy and Gaussianity and apply it to the Wilkin- 
son Microwave Anisotropy Probe (WMAP) cosmic microwave background (CMB) foreground reduced, temper- 
ature maps. We also test cross-channel difference maps to constrain levels of residual foregrounds contamination 
and systematical uncertainties. We divide the sky into regions of varying size and shape and measure the first 
four moments of the one-point distribution within these regions, and using their simulated spatial distributions 
we test the statistical isotropy and Gaussianity hypotheses. By randomly varying orientations of these regions, 
we sample the underlying CMB field in a new manner, that offers a richer exploration of the data content, and 
avoids possible biasing due to a single choice of sky division. In our analysis we account for all two-point corre- 
lations between different regions and also show the impact on the results when these correlations are neglected. 
The statistical significance is assessed via comparison with realistic Monte-Carlo simulations. 

We find the three-year WMAP maps to agree well with the isotropic, Gaussian random field simulations as 
probed by regions corresponding to the angular scales ranging from 6° to 30° at 68% confidence level. 

We report a strong, anomalous (99.8% CL) dipole "excess" in the V band of the three-year WMAP data and 
also in the V band of the WMAP five-year data (99.3% CL). 

Using our statistic, we notice the large scale hemispherical power asymmetry, and find that it is not highly 
statistically significant in the WMAP three-year data (< 97%) at scales ^ < 40. The significance is even 
smaller if multipoles up to ^ = 1024 are considered (~ 90% CL). We give constraints on the amplitude of the 
previously-proposed CMB dipole modulation field parameter. 

We find some hints of foreground contamination in the form of a locally strong, anomalous kurtosis-excess 
in the Q-i-V-i-W co-added map, which however is not significant globally. 

We easily detect the residual foregrounds in cross-band difference maps at rms level < 7/iK (at scales > 6°) 
and limit the systematical uncertainties to < 1.7/iK (at scales > 30°). 



I. INTRODUCTION 



Observational cosmology has established the flat ACDM 
model with nearly scale invariant initial density pertur- 
batio ns as the standa r d model o f mo d ern cosmology 
(e.g. lAstier et alJ (l2006h: ICole et alJ (l2QQ5l): lEisenstein et al 



(12005 ): iHinshaw et al.' ('2007'): 'Page et al.' ('2007'): i Riess et al 
(|2004i); [Spergel et al. (2007); Tegmarketal. (200i)). These 
observations seem consistent with the simplest predictions 
from inflation theory. Amongst those predictions, one con- 
sequence from the cosmological principle, the statistical 
isotropy (SI), and one generic consequence from inflation the- 
ories, the Gaussianity (to leading order) of the cosmic mi- 
crowave background (CMB) temperature fluctuations, have 
received a lot of attention with the release of the first year 
of observations of the WMAP satellite. The relevant sta- 
tistical analyses either aimed at detecting small amounts of 
non-Gaussianity (NG), that ste ms from non-li near effect even 
within inflation theories (Bart olo et al. 1 120041) . or looked for 
any anomalous signal that would challenge this standard 
model. 

However, separating SI from Gaussianity is a delicate 
task when making such a test, since one has to deal with 
only one realization of the CMB, that is considered in 
this context to be a random field. SI and NG have been 
tested in variety of ways and some "anomalies" have been 
reported. In particular, using tests opti mized for SI, in 
spher i cal ha rmonic (SH) phase space ( d e Qliveira-Costa et al.l 
1 19961 l2004l) an unusual alignment (98% CL) at low multi - 
poles have been found and confirmed (e.g lCopi et al.l (l2006h : 
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iLand & Magueijol (l2005ah ). Number of other tests and sta- 
tistical tools and estimators have been devised and used to 
constrain SI and/or NG. Among others, these i nclude: bi- 
polar power s pectrum (Hajian & Souradeepll2006h . phase cor- 
relations tests (iNaselsky et aLlbOOSh, higher order correlations 
in SH space ( bi/tri- s pectrum) e.g.('Cabena et al.| |2006l 120051: 
iFerreira et al.l 1 19981: iMagueii o & Medeiros 2 0041). n-point 
real s p ace statistics: dPurreret al.1 120001: iGaztanaga & Wag j 
I2OO3I: IGaztanaga et al.l l2003|), morphological estimators 
(like Min kowski functional) (IParkI [20041: IShandarinI 1 2002 : 
IWu etal." "200 1"), multipole vectors fCopi et al' l2006l 12004 : 
iLand & Magueijo 20 05b: Schwarz et al. 2004) higher order 
correlatio n functions (IGaztanaga & Wag3 2003|), ph ase space 
statistics (IChiang et al.l 120031: iNaselskv et al.l 12005 ). wavelet 
space statis t ics (ICruz et al.l 120071 : iMcEwen et al.1 l2 006a.b: 
Vielva et al. 120041) 7 ligher criticism sta tistic (ICavon et al] 
2005h . pair angular separation histograms (iBernu i et al.*' 2007h 
and also va r ious real- space base d tests eg: Eriksen et al.l 
(120071 l2004 : lHansen et aP (l2004allbh In particular a dedicated 
tests of hemispherical power asymmetry have been reported 
by many authors and found anomalous at confidence levels 
ranging from - 2cr to - 2.6cr (95%CL - 99% CL) 

In this work, we measure regional one-point statistics in the 
WMAP data and in simulations in order to test the SI and 
Gaussianity hypotheses. We mean to extend and generalize 
the previous similar works in three ways. 

Firstly, we show that the result of the analysis strongly de- 
pends on the way in which the sky is partitioned into regions 
for the subsequent statistics, and we circumvent this problem 
by relaxing the constraints on the shape and the orientation 
of a chosen sky pixelization by considering many randomly 
oriented sky regionalizations. This allows us to avoid a possi- 
ble bias in such regional analysis that is constrained only to a 
single choice of pixelization scheme. 

Secondly, we relax the constraint on the size of the re- 
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gions, thereby statistically probing features at different angu- 
lar scales. 

Thirdly, we account for all correlations between different 
regions, resulting from the well known two-point correlations 
(or possible higher-order correlations) using multivariate full 
covariance matrix calculus for more robust estimation of the 
statistical significance of local departures from Gaussian ran- 
dom field (GRF) simulations. 

We will assess the statistical significance of our results in 
three different manners so as to avoid the standard pitfalls of 
such an analysis and will rely heavily on realistic simulations 
to either probe the underlying distributions or to test the sen- 
sitivity of our statistic. 

The paper is organized as follows: in Sect. [Ill we introduce 
the data sets that are being tested, and provide details of the 
simulations. In Sect. Jill we describe the details of our statisti- 
cal approach for regional statistics. We then test and illustrate 
the sensitivity of our statistics via Gaussian and non-Gaussian 
simulations in Sect. HVlbefore presenting the results in Sect.lVl 
and discussing them in Sect. jVll We conclude in Sect. IVIII 



II. DATA AND SIMULATIONS 

For the main analysis in the paper we use the WMAP three- 
year foreground reduced temperature maps from differential 
assemblies (DA) Ql, Q2, VI, V2, and Wl, W2, W3, W4, pix- 
elized in the HEALPIX sphere pixelization scheme with res- 
olution parameter Ng = 512. We co-add them using inverse 
noise pixel weighting (Eq. [T]) and form either individual fre- 
quency combined maps (Q, V, W) or an overall combined map 
(Q+V+W) to increase the signal to noise ratio according to: 
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where Wi = X]j=j^i ^ji ^ji — ^ji/^oj the 
noise rms for a given DA and Nj j is the number of o bser- 
vations of the ith pixel for jth DA (iHinshaw et al.ll2007h . The 
sum over j iterates the DAs whose maps are co-added (in num- 
bers, {2, 2, 4, 8} respectively for Q, V, W and all channels). 
We will refer to those datasets as Q, V, W and INC (inverse 
noise co-added map) respectively and define a data set vector 
d G {Q, V, W, INC} for further reference. 

We also consider a difference maps between different chan- 
nels to independently test the residual foregrounds and to 
cross-check with the results obtained from the single band NG 
analysis. We consider a single band difference maps (e.g. Ql- 
Q2, VI -V2) as well, since nearly identical frequency differ- 
ence maps have a negligible amount of CMB or foreground 
signal^ and these are used to test the consistency of our white 
noise realizations against the pre- whitened 1 // pink noise of 
the WMAP data and constrain the systematical uncertainties. 



The non-vanishing CMB or foregrounds content, even in the single band 
differential maps, comes from slight differences in the effective working 
frequencies of the differential assemblies (DAs) and also from slightly dif- 
ferent beam profiles. While in case of the single band difference maps (e.g. 
Q1-Q2) the residual rms signal is weaker than the noise by more than two 
orders of magnitude, in case of the different frequency bands (e.g. Q-V) 
the residual CMB rms signal is about one order of magnitude weaker than 
the noise. 



Details of this check is given in appendix |Cl We will refer to 
these maps as QV, QW or VW for cross-band difference maps 
and Q12, V12 etc. for an individual differential assembly dif- 
ference maps. 

As an extension to the main analysis we also test the five- 
year WMAP data set from the V channel and refer to it as V5. 
For this purpose the WMAP five-year simulations are used and 
preprocessed in the same way as in case of the WMAP three- 
year data except for the sky-mask, which here we choose to be 
KQ75. 

The residual monopole, measured outside the three-year re- 
lease of the KpO (hereafter the Kp03 ) sky mask, is removed 
from each map by temperature shift in real space. The Kp03 
sky mask (including galactic region and bright point sources) 
is applied and no downgrading is performed at this level. We 
will use Nsim = 10^, realistic, full resolution simulations to 
test our statistics and to assess confidence thresholds (see Ap- 
pendix |Al for details and basic tests). 



III. DIRECTIONAL STATISTICS 

If the CMB sky is a realization of a multivariate Gaussian 
random field (GRF), then statistics of any linear statistical es- 
timator should not deviate from Gaussianity within any arbi- 
trary region in the sky. Otherwise - in case of non-linear es- 
timators - in general deviations from Gaussian statistics are 
expected, hence MC approach for assessing limits on consis- 
tency with Gaussianity is used. 

In order to test the stationarity and Gaussianity of the 
temperature fluctuation hypotheses we use two independent 
sphere pixelization schemes to define sky divisions and conse- 
quently a set of adjacent, continuous regions. 



A. Sky pixelizations 

The first pixelization scheme (hereafter referred to as HP) is 
an indep endent implementa tion of the HEALPIX pixelization 
scheme (iGorski et al.ll2005h and its resolution is parametrized 
by the Ug parameter (Fig. [T] top-left). The total number of 
pixels for a given is r = 12n^. We will use three different 
resolutions (n^) as specified in Table HI 

The second one (hereafter called LB) covers the sphere by 
dividing it along lines of parallels (iso-latitude) and meridians 
(iso-longitude) to obtain arbitrarily elongated pixels, generally 
of varying angular sizes (Fig. [T] top-right). This results in the 
total number of pixels r = NiNi) (where Ni and are the 
numbers of longitudinal and latitudinal divisions). The three 
different resolutions used in the analysis defined by these pa- 
rameters are specified in Table HI Further flexibility is allowed 
by rotating the polar axis by three randomly chosen Euler an- 
gles. 

Since there is no reasonable, physical motivation for pre- 
ferring any particular sky pixelization over another from the 
standpoint of testing a GRF hypothesis, we consider Nm = 
100 random orientations for each of the six types of pixeliza- 
tion schemes, which altogether yields 600 different sky pix- 
elizations with a total number of 280 000 regions of different 
shapes and sizes probing different angular scales (Fig.©. We 
therefore draw the three Euler angles used to define the axis 
position and pixelization scheme orientation about this axis 
from a uniform distribution. 
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Fig 1 : In the first row, two lowest-resolution pixelization schemes - 
HP 2 (top-left) and LB 32 8 (top-right) are shown. In the second row, 
we present an examples of two multi-masks actually used in analysis. 
Pixelization schemes are rotated to a random orientation, with Kp03 
sky mask applied. These are HP 2 (lower-left) and LB 64 8 (lower- 
right) respectively. Values in all regions were randomized for better 
visualization. 



All sky pixelizations are subject to the Kp03 (three-year 
KpO) galactic/point sources cut which masks ~ 23% of the 
sky. In practice there is no lower bound for the size of a region 
due to its random orientation with respect to the Kp03 sky cut. 
However for the sake of numerical stability, when computing 
the inverse covariance matrix (see below), we only consider 
regions that happen to have Npix > Npixth = 100, where 
Npix refers to the number of pixels of the original Ug = 512 
map falling into this particular region. 

Hereafter we refer to a particular random realization of a 
pixelization scheme (a random set of regions covering the full 
sky and merged with Kp03 sky mask) as a multi-mask, since it 
uniquely tags sky regions and allows to pursue statistics exclu- 
sively within them (see Fig. [T] bottom-left and bottom-right). 
Of course different multi-masks, even defined from a similar 
pixelization scheme, may have a different number of regions 
due to the random orientations with respect to the Kp03 sky 
mask. We define Nreg rn) as the number of regions of a 
multi-mask as a function of initial resolution parameter r and 
multi-mask ID number m G {l..Nm}- As an illustration, the 
two lowest resolution, pixelization schemes and two examples 
of multi-masks are shown in Fig.[T] We will also use additional 
sets of multi-masks to complete and extend the main part of the 
analysis in a few selected cases. 
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region No- 
Fig 2: Number of points, Npix{k), in regions for all 6 types of pix- 
elization schemes used in their initial position, after masking by 
Kp03 sky mask, in function of region number. All regions with 
Npix < 100 were not considered in the analysis as detailed in the 
text and treated as masked. The central parts (^ r/2) of a given 
pixelization are strongly covered by the Kp03 mask (this is true only 
for the particular initial orientation of a multi-mask). At the top ab- 
scissa we give the approximate angular scales probed by pixelization 
schemes with the corresponding total number of regions indicated in 
the bottom abscissa. 



skewness (S), and kurtosis (K)), of the underlying temperature 
fluctuations are computed for the data and for all A^sim = 10^ 
simulations. Together this yields 2.8x10^ regions for assess- 
ment of uncertainties. As will be shown in the Sect.|IVl allow- 
ing for arbitrary orientations of multi-masks has an impact on 
the results and yield a more stringent test on stationarity. The 
fact that we choose to work in real space allows for a good 
localization of deviations in the sky. 

The presence of extended, residual foregrounds or unre- 
moved, unresolved point sources, will affect the local central 
moments distributions. In particular the mean of the fluctu- 
ations will tend to be up- shifted with respect to simulations 
if diffuse foregrounds are present or down- shifted if they are 
over-subtracted. Also, depending on the amplitude of the 
residual foregrounds the local variance will also be altered. 
Looking jointly at the distribution of these moments on large 
scales might also provide a handle on the large scale distri- 
bution of power via the off-diagonal terms of the inverse co- 
variance matrix. The physical extent and position of the re- 
gions where particular type of deviation occurred can provide 
a clue to the possible nature of the foregrounds causing it (see 
SectlVl). 



B. One-point statistics 



In each of the defined regions of each multi-mask, the first 
four central moments (i.e. mean (m), standard deviation (a). 



Table I: Summary on the LB and HP pixelization schemes and resolu- 
tions used in the main analysis, given explicitly for quick reference. 
The columns abbreviations are as follows: (1) pixelization scheme 
reference name, (2) resolution parameter value, (3) approximated an- 
gular size of regions, (4) number of regions in pixelization scheme. 
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C. Assessing statistical significance 

Since our measurements are statistical, a crucial stage re- 
mains in probing their statistical significance. Our approach 
relies on a detailed comparison between the measurements 
performed on real data with the distribution of the same mea- 
surements performed on simulations. 

We consider three different ways to address the significance 
of these measurements. Each step involves one extra-level of 
generality and will shed light on the subtleties of such an as- 
sessment. 

At first, we look at individual regions, ignore their corre- 
lations and compare them with the simulations. We call this 
approach the "individual region analysis". It is the simplest 
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approach one can consider. 

Secondly, we compute the overall statistical significance per 
multi-mask, by taking into account the two-point correlations 
between moments of distributions (MODs) measured in re- 
gions of the same multi-mask via the full co variance matrix. 
We call this approach the "multi-region analysis". The result- 
ing probabilities ^(Xq) (Eq. IB7I) are the joint probabilities 
of exceeding a certain confidence threshold as a function of 
pixelization scheme (r), multi-mask (m G {l..Nm}), MOD 
(X G {m, cr, S, K}), and dataset (d G {Q, V, W, INC}) con- 
figured by a parameter vector q = {X, r, m, d}. This anal- 
ysis extends the information from the single region analysis 
by testing the consistency of the data with the simulations via 
standard multivariate calculus. 

Finally, we combine all the information probed by different 
multi-masks to find the joint cumulative probability of reject- 
ing the GRF hypothesis as a function of pixelization scheme 
(r), MOD (X), and dataset (eg. frequency) (d). We call this 
approach the "all multi-masks analysis". 

We remind that the statistical significance of any real data 
measurement, at any stage of the analysis, is always assessed 
by a comparison to the set of the same measurements per- 
formed using GRF simulations. The exact details of the anal- 
ysis at each step are given in Appendix |Bl 

D. Visualizing the results 

To visualize our results from the single-region analysis, or 
multi-region analysis at certain confidence level, we proceed 
the following way. For individual region statistics, for each 
region of each multi-mask we define Uo- as 

= V2eTr\l - P{X)) = cdfG"'(P(X)/2) (2) 

where P{X) is the quantile probability derived according to 
Eqs. lB3l and lB4[ The Ua- thus defined is the Gaussian number 
of as by which a region, defined by a given multi-mask, de- 
viates from simulation average. We then produce maps of Ua- 
estimator, for data processed through each of the 600 gener- 
ated multi-masks and for each MOD. Then, to present all the 
results in a compact way, we scramble these maps within the 
same MOD. We over-plot the individual pixels from regions 
with the strongest deviations from the underlying pixels. Pos- 
itive Ua- values correspond to excessive value of a given MOD 
in a region, and negative value correspond to its suppression. 
For clarity, we use a threshold \ncj^th\ = 3 to produce maps 
with only the strongest (3a) detections. 

For the joint multi-region statistics we produce maps 
(as detailed above) using only those multi-masks that yield 
|P(Xq)l ^ Pth (Eqs. lB3llB4l) revealing detections at the sta- 
tistical significance 1 — Pth for a given MOD, multi-mask res- 
olution r and dataset d, i.e. for a given parameter q. Within 
this notation, single region statistics corresponds to Pth = 1. 

While the "n-sigma" maps are easy to read when looking 
at distributions of the anomalies at a given local significance, 
they are unitless and cannot be directly linked with quanti- 
ties that are physically measured. We therefore also consider 
a difference maps (A maps^) of regional departures in indi- 



^ We will hereafter refer to maps so produced as A maps to make a clear dis- 
tinction from the difference maps obtained by differentiation of temperature 
maps from different frequencies (e.g. QV, QW, etc..) 



vidual MODs between datasets and averaged simulation ex- 
pectation: i.e. for ith region of a given multi-mask we plot 
Ai = Xi — (Xi) AT, where () at stands for average over N sim- 
ulations. 



IV. TESTS OF THE SIMULATION AND MEASUREMENT 
PIPELINE 

In order to validate the statistical tools introduced above, 
test the sensitivity and the correctness of the numerical code, 
we performed a set of experiments using both simulated 
WMAP CMB data and data with either simulated violation 
of the large scale statistical isotropy or localized NG features. 



A. Consistency check with GRF simulations 

We first test self-consistency by generating 10 additional 
INC CMB data sets and carry out for each of them the single- 
region, joint multi-region and all multi-masks statistics. As 
expected, we found that the simulated datasets yield a good 
consistency with the simulations (at 68% CL). 



B. Sensitivity to local NG 

For a statistically isotropic Gaussian process the kurtosis is 
expected to be exactly equal to K = 3, which translates into a 
kurtosis excess KE = K — 3 = 0. Violation of either of the 
assumptions can lead to a positive or negative KE. 

At first we simulate what could be a residual component, re- 
sulting from subtraction of an non-ideal foreground template, 
extended over an area of 10° angular, centered at = 
(50°, 50°), whose signature would be a non- vanishing KE. 
Such residuals must be expected to be small in the foreground- 
cleaned maps, and they could be either positive - resulting 
from foregrounds undersubtraction - or negative ones - re- 
sulting from foregrounds over subtraction. 

Therefore, the introduced NG component is drawn from a 
normal distribution with variance and with a mean increas- 
ing uniformly across the patch when going north-south follow- 
ing Healpix ring ordering inside the spot. From the northern 
point of the patch down to the southern part of the patch, the 
mean will shift by 2nacMB- 

We introduce in this way a gradient in the noise and thus a 
non-zero local negative KE (Fig. [4]) but still preserve a vanish- 
ing (within the spot) skewness. The a value is chosen to be 
1% of the underlying CMB rms (ctcmb). To test the sensitiv- 
ity of our estimator we consider n to be either 1 or 2, above 
which the NG template starts to be visually noticeable due to 
edge discontinuities. Note that, with a so defined anomaly, the 
pixels of the spot that are close to its horizontal diameter will 
have the least impact on the underlying CMB field distortions. 

The choice of the n parameter values corresponds to the 
NG signals of the rms amplitude ~ 50/iK and ~ 100/iK for 
n = 1 and 2 respectively within the spot. We note that the 
rms. value within the spot of the same size, in the foregrounds 
reduced difference VW map of the WMAP3 data, yields > 
100/iK depending on the exact location of the spot in the sky, 
hence our choice of NG signals amplitude aim at detection of 
relatively small anomalies as compared to the WMAP3 noise 
specifications. 
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Fig 3: aj A Non-Gaussian temperature gradient template, leading to 
locally negative kurtosis excess. The parameter value of n = 2 is 
used (see the main text for details). Note that one of the spots is 
practically removed by the galactic sky cut. b-d) The results of the 
single-region analysis using the 100 multi-masks of the HP 4 pix- 
elization scheme of one of the simulated INC data and contaminated 
with the template. The scrambled ria detection map in variance, 
skewness and kurtosis, thresholded at na,th — 3, is plotted. Note 
the strong, template-induced, local anomalies detections traced by 
different multi-masks, as well as some other, but somewhat weaker, 
detections of | rial > 3 regions. In particular, the template leads to 
the local kurtosis suppression and strong local excess of the variance. 



We find that for a single NG spot of radius 10°, the multi- 
region analysis does not return any significant detection for 
n = 1 in any of the MODs, but the single region analysis 
finds the contaminated regions unusual at ^ 2.5. In case 
of n = 2 we detected a 3cr local deviation in kurtosis, while 
in ail-multi-masks analysis we reject Gaussianity at 99.8% CL 
(HP 8) due to variance distributions^. 

We rerun the test for the same type of NG templates but 
replicated in 3 disjoint spots at different directions in the sky 
for the same values of n parameter (Fig.[3]-a). 

The choice of the position and the size of the spots is most 
relevant to the results presented in Section |Vl The results of 
the single region analysis is shown in Fig.[3b)-(d). Note how 
different multi-masks trace the locally introduced anomaly. 
Depending on the orientation of the multi-mask and its regions 
around the directions of the NG spots, the returned Ua- val- 
ues differ. In the overall multi-region analysis this naturally 
leads to a distribution of probabilities which strongly depend 
on how the features of the map are split and probed by differ- 
ent regions. 

Note that some of the multi-masks also return an > 3 
detections even in a template-free regions. It is therefore clear 
that use of many differently oriented multi-masks helps to in- 
vestigate the statistical significance of local anomalies. 

The multi-region analysis in case of n = 1 return no sig- 
nificant detections in any of the MODs, but very significant 
deviations were detected for the case n = 2 in dXl-multi-masks 
analysis (Table|IIlin section "KE-"), again only in the variance 
distributions. 

Consequently, we find that for the unfiltered maps the dis- 
tributions of variances are actually more sensitive to this kind 
of simulated anomalies, rather than higher order MODs. We 



note, however, that measuring a local sign of KE may be a hint 
of the nature of the foregrounds signals as the kind of template 
used in this example introduces locally only the negative KE 
as shown in Fig.lH Similar dependences are obtained for tem- 
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Fig 4: Negative departure of the kurtosis in the NG spot as a function 
of n parameter. The plot shows the kurtosis of a sum of GRF and 
the generated NG template averaged over 10 random realizations of 
the two. The variance of the GRF is 100 larger than the variance 
of the NG template and the field size is chosen to correspond to the 
NG spot size as described in the text. The 3cr confidence thresholds 
were quoted around the expected values for the GRF only and for the 
GRF+NG template cases. Vertical lines indicate the chosen values of 
the n parameter used for the templated maps generation for the test. 

plates of different shapes and sizes and combinations of these 
factors. In the limit of a flat template (n = 0) the field be- 
comes Gaussian, as expected. 

It is possible to introduce a non-vanishing skewness by 
making the template unsymmetrical, by considering shifts in 
the mean which are unsymmetrical about zero. A similar 
effect is possible by considering regions of multi-masks that 
only partially overlap with the area of the NG spot (Fig. [3]:). 
Note that the example from Fig. (H does not include the ef- 
fects of the non-uniform noise component which is included 
in our tests and also that the confidence level contours were de- 
rived assuming the Gaussian error statistics; therefore it is not 
straightforward to extrapolate the strong detections expected 
from Fig. (H for n = 2 onto the full sky, locally templated, 
signal and noise realizations, subject to subsequent regional 
statistics. 

Note that the assumed size of local deviations is small rel- 
ative to the full sky one, and hence their global impact is re- 
duced accordingly. Larger, in sense of area, deviations will be 
of course easier to detect. Also a specific pre-filtering in the 
spherical harmonic (SH) space, of the data prior to the test may 
help to expose the most relevant scales to the test. In this test 
we focused on testing unfiltered maps; therefore, necessarily 
the strength of the detections must be suppressed. 

It is not possible to obtain a locally positive KE with the 
above-described template. Such deviation could however be 
the signature of the unresolved point source contribution"^, or 
an unknown and localized noise contribution. 

Again, for qualitative studies only, we simulate the point 
source component in the full sky by adding random num- 



^ Of course, the estimated rejection confidence thresholds given in Table. |lTl 
based only on a single- simulation measurements may be somewhat biased 
depending on particular realization of the GRF simulation. 



Although this would have a specific frequency dependence, we ignore this 
fact here. 
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bers drawn from a distribution which is an absolute value of 
the Gaussian distribution with zero mean and random vari- 
ance nacMB parametrized by parameter n and uniformly dis- 
tributed within range [0,ncrcMB], where acuB is the rms 
value of the underlying CMB fluctuations. 

From our simulations however, it appears that it is difficult 
to detect a significant contribution due to point source contam- 
ination since such signal is significantly smeared by the instru- 
mental beam even for n as large as 6. Even when KE > 6 
before beam smearing, the variance response is much stronger 
than the KE response leading to inconsistencies with simula- 
tions in the total power of the map as measured by e.g. full 
sky variance distribution (Fig. [18]). We therefore conclude that 
it is unlikely to detect any point source contribution to KE in 
this test which is not surprising since we work at fairly low 
resolution which dilutes the point source signal. 

As already mentioned, a locally generated noise-like com- 
ponent in the map, in principle, could generate an non- 
negligible positive KE^ as it would not be processed by the 
instrumental beams. However since such noise is not well mo- 
tivated physically due to non-local properties of the TOD data 
and scanning strategy of the WMAP, and also since the noise 
properties are well constrained, therefore we do not consider 
such case. 

We conclude that small and single (compared to full sky 
observations) localized NG features will be difficult to detect 
via higher order MODs in the joint multi-region and dXl-multi- 
masks analysis due to their small statistical impact on the over- 
all statistics. However a single region statistics carried out first 
might be a rough guide in selecting a possibly interesting fore- 
ground NG signals. If these indicate regions with significant 
deviations in variance and possessing a negative KE it would 
hint on residual, large scale foreground contamination. 



1. Stability of results in function number of multi-masks 

As different multi-masks probe the underlying data differ- 
ently, the joint-probabilities per multi-mask differ and lead to 
a distribution that typically covers a wide range of possible 
probability values. As such, the multi-region analysis (see 
sec. IB 21 for details) can be used to find the orientation of the 
most unusual regions in the multi-masks that yield the smallest 
probability as comparred to GRF simulations. 

Since from the point of view of statistical isotropy all multi- 
masks are equivalent, in the dX\-multi-masks analysis (see 
sec. IB 3l for details) we integrate the infomation from all multi- 
masks within a pixelization scheme to obtain an average level 
of consistency of the data with GRF simulations for that pix- 
elization scheme. This approach also provides a conserva- 
tive way of averaging over a possibly-spurious detections that 
could be a fluke, due to some accidental arrengement between 
a multi-masks and a data set. If the anomalous feature in the 
map is strong enough to be detected in many multi-masks then 
this will also result in a detection in the joint dXl-multi-masks 
analysis. Conversly, if only one or few multi-masks result in a 
very small probability the overall impact will not be large due 
to stability of the median estimator with respect to the distri- 
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^ This is most easily seen in the contribution of the anisotropic noise of the 
WMAP to the kurtosis of the signal only simulated map, which induces a 
significant positive overall KE value in the simulation. 



Fig 5: Convergence of the joint all multi-masks probability to the 
value reported in Table |TT1 (in section (2KE-) for n = 2) for the case 
of HP 2 pixelization scheme in function of the cumulative number 
of multi-masks used to derive it (black solid lines) for mean (top 
panel), variance, skewness and kurtosis (bottom panel). Addition- 
ally we overplot the joint multi-region probabilities per multi-mask 
for each multi-mask number (blue dashed lines). 



bution outliers. However, the investigation of the most anoma- 
lous multi-masks may help in selecting the deviating regions 
for further analyses. 

In this section we show how the convergence to the results 
of the Sill-multi-masks analysis is reached in function of num- 
ber of multi-masks used to derive the median value and 
the corresponding median -distribution leading to the all- 
multi-masks probability. As shown in Fig. [5] the convergence 
of the joint dXl-multi-masks probability to the reported value 
for = 100, is in this case {HP 2 pixelization scheme) quite 
fast; however in general it depends on particular properties of 
the map as well as on the set and ordering of the multi-masks 
used. The speed of the convergence and the robustness of the 
final value is as good as the convergence of the unbiased mean 
estimator - i.e. the median - to the intrinsic mean value as 
the number of random trials (corresponding to the number of 
multi-masks) probing the underlying distribution increases. 

Note how strongly the joint probability per multi-mask de- 
pends on the orientation of the multi-mask (blue-dashed lines 
in Fig. [5]). In particular for the variance case, the probabilities 
of rejection per multi-mask range from 45.7% for multi-mask 
m = 51, to 99.8% for multi-mask m = 2. Yet the, reported 
in Sill-multi-masks analysis, median value is very stable with 
respect to these variations. 



2. The shape of multi-mask regions. 

Naturally, the multi-masks having large regions will be in- 
sensitive to the small scale map features, while the multi- 
masks having small regions will be insensitive to the large 
scale structures. This motivates the usage different number 
of regions to probe different scales of the map. 

Theoretically there are infinitely many ways of defining the 
shape of regions of multi-masks, and of course, our choice 
of the shape of pixelization schemes and the multi-masks is 
somewhat arbitrary; however the motivation for using differ- 
ent shapes of regions is straightforward: to probe the data 
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using different binning techniques since the GRF statistically 
should not depend on it. However if the data turns out to be 
non-random it is possible that the non-randomness will be ex- 
plored differently by different region shapes. A loose analogy 
to the real-space multi-mask region shape, which is used to 
derive a local value of an estimator, over the defined area for 
a given orientation of multi-mask, is in the wavelet space the 
shape of the mother wavelet, which is used to obtain the local 
convolution coefficients for the input map. In this analogy the 
size of the region correspons to the wavelet scale. 



C. Sensitivity to the large scale phase anomalies 

We test the sensitivity of our method to the well known large 
scale anomalies found in the WMAP data: i.e to the aligned 
and planar low multipoles £ = 2 and £ = 3. In order to 
test such anomalies we generate two GRF CMB simulations. 
In both simulations we use the same realization of the power 
spectrum and phases as in the case of the first GRF simulation 
rSect. HVAl) . 

In the first simulation we enforce large scale phase corre- 
lation by introducing an "m-preference" in the power distri- 
bution in the quadrupole (£ = 2) and octupole (£ = 3). We 
choose the "sectoral" spherical harmonic coefficient an, to 
carry all the power of the multipoles. In the second simula- 
tion we extend this modification up to ^ = 5. The GRF signal 
simulations are rotated to a preferred frame prior the introduc- 
tion of the planar multipoles. The signal maps are then rotated 
back to the original orientation so that the maximal momen- 
tum axis was located at = (260°, 60°) before adding 
noise. 

As shown in Table HIl such anomalies have not been signif- 
icantly detected at tested scales. Although it is expected and 
observed that considering unfiltered maps (containing all mul- 
tipoles information mixed together) there will be a little over- 
all impact on the statistics we note a higher sensitivity would 
be obtained is a prefiltered in SH space data were used. 



Table II: PA\-multi-masks analysis results from tests of the simulated 
datasets. The columns content is as follows: (1) pixelization scheme, 
(2.. 4) "probability of rejecting" the consistency with GRF simula- 
tions fSect. lB 3l ) for each MOD. The abbreviations of the datasets are: 
1KE-: INC simulation with one NG spot (three NG spots) leading to 
KE< and n = 1 (see text Sect.|TVlfor details), 2KE-: INC simula- 
tion with one NG spot (three NG spots) leading to KE< and n — 2 
(see text Sect. |lVl for details), A2..3: INC simulation with aligned 
multipoles £ — 2 and £ — 3, A2..5: INC simulation with aligned 
multipoles from £ = 2 to £ = 5,M: INC dipole modulated simula- 
tion with dipole amplitude of Ai^^^. <M>: V simulation with CMB 
signal fully (partially) modulated according to parameter Ai^^^ . The 
average confidence thresholds are given. The values are rounded to 
integer percentiles in case of probabilities < 99% CL. The saturated 
values are marked with 

Rig^JE p;2.[%] p^~m 

(1) (2) (3) (4) (4) 

n=l: IKE- (1KE-): 
no significant detections (no significant detections) 



n = 2: 2KE- (2KE-): 



HP! 


42 (45) 


82 (94) 


58 (75) 


49 (80) 


HP A 


43 (69) 


97 99.98) 


31 (66) 


39 (62) 


HP^ 


8(70) 


99.8 (99.8) 


5(32) 


23 (41) 


LB 32 


28 (58) 


88 (99.0) 


29 (72) 


27 (53) 


L5 64 8 


34 (64) 


69 (96) 


20 (60) 


18 (34) 


LB 64 16 


18 (69) 


98 (99.5) 


10 (37) 


28 (44) 






A2..3 (A2..5): 








no si| 


^nificant detections (no significant detections 




M: Aio24 = 0.114 (Aio24 = 0.2) 




HP 2 


39 (36) 


94 (99.99) 


42 (50) 


32 (35) 


HP 4 


43 (47) 


99.6 99.99*) 


22 (25) 


31 (32) 


HP^ 


11(19) 


99.96 99.99*) 7(8) 


18 (20) 


LB 32 S 


30 (35) 


97 99.99*) 


20 (22) 


22 (23) 


LB 64 ^ 


34 (38) 


91 99.99*) 


15(17) 


14 (14) 


LB 64 16 


22 (30) 


99.4 99.99*) 


12(13) 


23 (25) 




<M>: Aio24 = 0.114 (A40 


= 0.114) 




HP 2 


58 (59) 


99 (56) 


52 (44) 


54 (44) 


HP 4 


77 (56) 


99.87 (57) 


54 (37) 


50 (38) 


HPd, 


79 (62) 


99.97 (55) 


64 (29) 


54 (52) 


L5 32 8 


78 (57) 


99 (55) 


54 (38) 


51 (47) 


LB 64 ^ 


78 (58) 


97 (55) 


59 (38) 


54 (43) 


LB 64 16 


76 (61) 


99.7 (54) 


64 (32) 


55 (46) 



D. Sensitivity to the large scale power anomalies 

We give a special attention to testing the sensitivity of 
the method for detecting and quantifying the previously re- 
ported large sc ale an omaly in the power distribution in the sky 
dEriksen et all2007ll2004i) . We create a simulated CMB maps 
where the CMB signal is modulated according to: 

T(n) = TcMBin){l + M{n)) 

M(n) = A,_n-d ^ ' 

where fi is a unit vector and M is a bipolar modulation 
field, oriented in direction d = (225°, —27°) with amplitude 
^^max ^ {0.114, 0.2} which modulates the CMB component 
up to the maximal multipole of ^max = 1024. 

The result of the test with such modulated simulation is 
given in Table HIl (in section "M"). As the amplitude of A = 
0.114 has been previously cl aimed to be preferred for the 
CMB data dEriksen et al.ll2007h we process five additional, full 
resolution V simulations, modulated with this amplitude, to 
reduce a potential biases from a single draw of a random sim- 
ulation, and report (Table HIl in section "<M>") the average 
rejection thresholds as a function of the pixelization scheme. 



We also process five additional modulated simulations, and 
apply the modulation only to the range of multipoles ^ < 40, 
leaving higher multipoles unmodulated, since the aforemen- 
tioned work operated at much lower resolution. 

The test is able to reject the modulation of the CMB of 
the amplitude A1024 = 0.114 at a very high confidence level 
(99.9%) depending on the pixelization scheme. Note however, 
that this model modulates all scales equally. 

Although in principle, the modulation will change the un- 
derlying power spectrum at scales where it was applied, we 
estimate (Appendix IA 1 al) that any such effect for the modula- 
tion of A1024 = 0.114 still remains in greatly consistent with 
the non-modulated simulations' power spectrum, and hence 
the results given in Table HIl do not result from the underlying 
power spectrum discrepancies. 

We are unable to reject the possibility of modulation with 
such amplitude applied only to the large scales (^ < 40). Such 
modulation is therefore consistent with the GRF field or unno- 
ticed by the test (Table HJ). However, according to the best-fit 
ACDM model (assuming even a noiseless observation), the 
multipoles £ < 40 carry only about 24% of the map's power. 
The possible modulation signals at these scales must also be 
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Fig 6: Histogram of the reconstructed modulation axis orientations 
as measured via the minimal joint probability P{x^) fEq. IB7I) in 
the set of 96 LB \ 2 multi-masks from 1000 modulated simulations 
(top-left) visualized using Healpix grid of resolution Us = 4. The 
corresponding, reconstructed 50% (dark blue), 68% (light blue) and 
95% (light red) confidence level contours obtained after smoothing 
the histogram with a Gaussian beam of FWHM= 7° (top-right). The 
reconstructed distribution is normalized on a hemisphere. In the bot- 
tom plot: the distribution of minimal log-probabilities (log P (xq) ) 
obtained from 1000 simulations modulated with amplitude A40 = 
0.114 along with the value obtained from the V map of the WMAP3 
data. 



more difficult to constrain as these are dominated by the cos- 
mic variance uncertainty. 

To investigate this further, we test simulations with only the 
large scales being modulated according to A40 = 0.114 along 
direction = (225°, -27°). We filter out these scales 

up to fmax= 40 (using the Kp03 sky cut) in SH space and 
downgrade the map to rig = 64, and process these using a 
new set of multi-masks of type: LB 1 2 - i.e. having only two 
regions, each covering a hemisphere. We use 96 such multi- 
masks with orientations defined by the centres of pixels of the 
northern hemisphere in the ring ordering of the Healpix pix- 
elization scheme of resolution Us = 4. We prepare a set of 
1 000 modulated simulations treated as data and use 1 000 in- 
dependent GRF simulations to test the consistency with SI. 
We split the GRF simulations into two sets of 500 simulations 
each, to derive the covariance matrix, and probe the under- 
lying distribution.^ We carry out the multi-region analysis 
using 96 multi-masks LB 1 2, and record the values of minimal 
probabilities (per modulated simulation) and the correspond- 
ing orientation of the multi-mask. The spatial distribution of 
these orientations defines the accuracy to reconstruct the cor- 
rect intrinsic modulation field orientation at these scales, under 
cut sky and negligible amounts (at these scales) of noise. 

The result is plotted in Fig. [6] (top-left). It is easily seen that 
while the direction is quite correctly found, the dispersion of 
the directions within even 50% CL contour (Fig. [6] top-right) 
is quite large which precludes a very precise determination 
of the modulation axis. We find, that statistically ~ 8% (the 
probability corresponding to the peak value in the bottom plot 
in Fig. [6]) of Gaussian simulations, to which we compare the 
modulated simulations treated as a data, exhibit a more un- 
usual configurations of hemispherical variance distribution. 



^ Note that with the LB I 2 multi-masks having only two regions it is not 
necessary to process a very large number of simulations to assess a good 
convergence. 



mean variance 




-3.7 0.1 3.afn^l ^4.? -0-0 :U{tK,] 



skewness kurtosis 




^.S 0.3 riA\n^] '4.4 0.^ ,'i.2|i(^J 



Fig 7: Results of the single-region analysis visualized in the compos- 
ite ria maps of the INC data, for each of the MODs. The threshold 
of na,th = 3 (as defined in Sect. lIIIDl is used. Mollweide projection 
and the galactic coordinates are used. The origin of the coordinate 
system (/, 6) = (0, 0) is in center of the plots and the galactic lon- 
gitude increases leftwards. Regions around the Galactic plane are 
partially removed by the Kp03 sky mask. The same convention is 
kept throughout the rest of this paper. 



Consequently, we report this median rejection threshold of 
~ 92%, as the statistical sensitivity of the method for detect- 
ing the large scale modulation (i.e. in the filtered maps mod- 
ulated with A40 = 0.114). We consider this result - i.e. the 
low rejection confidence level - to be penalized mostly by the 
cosmic variance uncertainty and freedom of phases to assume 
an arbitrary orientations with respect to the (unsymmetrical) 
sky cut. Consequently, we note that due to these uncertainties, 
it may be difficult to increase this rejection level for scales of 
^max~ 40 and amplitude A = 0.114. 



V. APPLICATION TO WMAP THREE-YEAR DATA 

We now present the results of the statistics described in 
Sect. Unl applied to the three-year and five-year WMAP data. 



A. Individual region statistics 

The individual region statistics as described in Sect. Hill find 
numerous regions amongst our many pixelization schemes, 
which deviate by more than certain \na^th\ in all MODs. Ta- 
ble |llll gives an incomplete list of some of the strongest detec- 
tions. In Fig. [7] we plot the detected regions at {ria^thl = 3 in 
the individual region statistics of the INC data. 

Note that with our approach the so-called "cold spot" is not 
detected at 3a level from being cold (i.e. via distribution of 
means) at all tested resolutions (Table Hill). It appears at about 
~ 2.7a around galactic coordinates = (211°,— 57°). 

Excessively "cold" or "hot" deviations in all MODs are de- 
tected in general with large values of Ua- . 

As expected the Ua- map of the means in Fig. |71 shows that 
the strongest deviating regions are directly close to the galactic 
plane cut off by the Kp03 mask, thus hinting at foreground 
residuals. 

The variance Ua- map shows local strong anomalies with 
the extended variance suppression in the northern hemisphere 
towards (/, b) = (67°, 19°) and with an extended variance ex- 
cess towards (/, b) = (199°, -55°). We note that these local- 
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Fig 8: Results of the single region analysis. Residual difference 
(A = \fo^ — (a/o^)) scrambled map of the variance distribution 
in the INC map processed with 100 multi-masks of the HP 2 pix- 
elization scheme. The well know, large scale hemispherical power 
distribution asymmetry is clearly seen as is the distribution of fore- 
ground residuals along the galactic cut. 



ized anomalies must, at least in a part, make up for the hemi- 
spherical power asymmetry. 

Skewness and kurtosis maps consistently indicate strong 
local deviations from GRF simulations towards (/, 6) = 
(193°, -26°) and (/, h) = (356°, -36°). While some regions 
appear in all three maps, some appear only in one of the mo- 
ments therefore the correlation between those results is not 
obvious. 

In order to investigate the spatial distribution of this high 
Ha- regions, we plot in Fig. [8] the A map, of differences be- 
tween the variance distribution measured in each region of 
multi-masks of the HP 2 pixelization scheme, and the sim- 
ulations average, A = di = — {^/~a^). The map is 
obtained by scrambling 100 difference maps from 100 differ- 
ent multi-masks (as described in Sect. lIIIE)!) . This map can be 
seen as a residual map for the local variance. This residual 
map exhibits a well known power asymmetry dEriksen et al.l 
[2007, 2004j). The dipole (^ = 1) component of this residual 
masked map (ignoring the effect of the mask) is aligned along 
axis (/, 6) = (237°, —44°) with power excess in the southern 
hemisphere. In order to probe this direction further, for each 
individual multi-mask we also produced an n^j estimator maps 
and checked the orientation of the dipole axis. Fig. [9] shows 
the PDFs obtained for the orientation of the dipole asymmetry 
in galactic longitude and latitude. Interestingly, we notice that 
the orientation of the axis of the hemispherical power asym- 
metry has some scale dependence. When using smaller scales 



Table III: An incomplete list of strongest deviation directions from 
maps in Fig [7] at \na,th\ = 3 in the individual region analysis for 
means, variances, S and K, sorted in galactic longitude ascending 
order. Notice that with our simulations we directly probe the ±3.7cr 
PDF region without need for extrapolations. 

mean variance skewness kurtosis 

(l,b) rig (l,b) ricj (l,b) riq (l,b) 

157,-29 -3.07 63,28 -4.21 173,-73 -4.50 84,-31 -4.38 

211, -57 -2.78 67, 19 -4.21 193, -26 5.06 195, -27 5.20 

241,42 -2.36 199,-55 4.13 209,8 -3.50 212,-55 -3.77 

265,-21 3.81 319,-27 3.40 217,35 4.00 309,59 -3.25 

318, -9 -3.69 225, -20 -3.60 312, -21 3.65 

167, 79 3.06 241, -52 -3.51 356, -36 4.25 

311,-20 3.86 

357,-35 4.15 




180 190 200 210 220 230 240 -70 -60 -50 -40 -30 -20 -10 

I [deg] b [deg] 

Fig 9: Orientation of the dipole component (galactic longitude - left 
panel and galactic latitude - right panel) of the ria maps for the INC 
map processed individually with all HP pixelization schemes. Each 
color corresponds to a different resolution of the pixelization scheme. 
While the longitudinal orientation of the dipole does not vary with 
resolution, the galactic latitude systematically shifts to lower galactic 
latitudes as resolution increases. 



with finer pixelization schemes, the orientation of the power 
asymmetry dipole shifts from larger galactic latitudes (roughly 
from the position of the cold spot, with the mean PDF value 
(/, h) = (218°, -43°) - see also Table|IIIl) for the HP 2 reso- 
lution to smaller latitudes (/, h) = (206°, -18°) for the HP 8 
pixelization scheme. The dependence of the power asymme- 
try orientation in function of th e pre-filtered i n SH sp ace data 
have previously been tested by iHansen et al.l (l2Q04ah . While 
we will return to the power asymmetry issue in the next sub- 
sections, we note that the medians of the dipole axis distribu- 
tions of other MODs are not correlated with the dipole axis 
orientation of the variance map (Fig.O and generally point at 
some other locations. 

In the next section we quantify the statistical significance of 
these deviations. 



B. Joint multi-region statistics 

In Fig.[TOlwe plot a compilation of all joint "probabilities 
of exceeding", calculated with all datasets considered (Q, V, 
W, and INC) in all 600 multi-masks. In order to visualize 
the smallest probabilities logarithmic scale is used. Note that 
we sort these probabilities in each MOD and dataset so as to 
ease visualization, so the points with same abscissa in differ- 
ent MOD and data sets do not necessarily correspond to the 
same multi-mask. 

Most of the results concentrate along the zero point of the 
joint log-probabilities, which indicates a good consistency of 
the data with the simulations at relatively high CL. (The white 
region in the Fig. [TOl encompasses CLs of up to 3cr; 3a and 4cr 
regions are shaded in red and yellow respectively). 

It is important to note that within one pixelization scheme 
the dispersion of probabilities in the Fig. [TOl results only 
from the orientation of the multi-mask. As a result, the 
statistical method involving many pixelization schemes help 
us obtain the unbiased results that one could get relying 
only using a single pixelization scheme. We also recall that 
for each plotted point, the statistic was also calculated for 
A^simPDF(x2 ) = 1 000 simulations in order to probe the 
underlying PDFs. For each point the corresponding full 
covariance matrix was obtained from A^sim(Cq) = 9 000 
simulations (as described in Appendix IB 21) . 

We now focus on three distinctive sets of results, based on 
the Fig. [TOl and quantify the deviations in more detail. We 
detail on a tentative excess seen in kurtosis before focusing 
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Fig 10: Results of the "multi-region" analysis. Dependence of joint probabilities of exceeding (P(xq), Eq. lBVI l as a function of multi-mask 
number for all pixelization schemes considered. The probabilities calculated with the spectral data maps from channels Q, V, W, INC are 
plotted in red, green, blue and black dots respectively. Each dot corresponds to the joint probability using one multi-mask. From the left to the 
right, the panels show results for increasing resolution of the pixelization scheme (see Table U) with 100 different multi-masks (along abscissa) 
in each. From top to bottom the four rows correspond to the four MODs - i.e. mean, variance, S and K respectively. Probabilities corresponding 
to the WMAP5 V5 data for pixelization scheme HP 4 are plotted with green crosses (+). 3a and 4a confidence levels are shaded in red and 
yellow respectively. The joint probabilities were sorted in each dataset before plotting for better visualization; therefore the probabilities from 
different datasets generally do not correspond to the same multi-mask numbers and do not directly correspond to the unique reference numbers 
used in the analysis. Hence the most unlikely events are localized at the left side in each panel. Additionally we plot the thin red line, which 
indicates the distribution of probabilities obtained from 100 GRF simulations of the Q data, each of which was processed with one (different 
for each simulation) multi-mask. If the data follow the expectation of GRF then statistical departure of data from this line would manifest 
certain degree of correlation between probabilities obtained with different multi-masks for the same dataset (as discussed in Sect. rvB . For 
better visualization in range m G [10, 90] we plot only every lO'th sorted probability value. 



on the large scale power anomalies, and on unusually strong 
dipole contribution in the V channel of the WMAP. Then, we 
comment on the results from tests carried out with the differ- 
ence map datasets. 

1. Localized Kurtosis excess 

In Fig. [TOl there is one "3cr" detection in kurtosis in HP 4 
pixelization scheme in the INC dataset (bottom second from 
the left panel in Fig. [T0|) - a result found using one in 100 of 
multi-masks probing these scales. Here we discuss this partic- 
ular point as a tentative detection because although the multi- 
mask bins the data to create the most unlikely realization of 
the kurtosis, it lies in the low-end tail of the whole spectrum of 
equivalent measurements and hence its statistical impact can- 



not be large. 

Table IV: Three a NG detections in K multi-region statistics of the 
INC data using multi-mask resulting in joint probability P{Xq) < 
0.0027 for the resolution Us = 512 and P{xl) < 0.005 for the 
resolution ris = 64. The (l,b) field gives the galactic coordinates to 
the center of the region 





resolution (us 


= 512) 


resolution (ris 


= 64) 




region 




Ua 






(lb) 


160 


2750 


3.48 


37 


1.38 


181 ,2 


185 


12610 


3.21 


199 


3.26 


199 , -23 


104 


16125 


3.52 


250 
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In Fig.fTTh we plot, thresholded at "3cr", the Ucr map using 
only this particular multi-mask. The details of the three most 
strongly deviating regions in this map are given in Table [IVl 
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Fig 11: a) Kurtosis (thresholded at 3a) estimator map from the 
multi-region analysis of the INC data, inconsistent with GRF hypoth- 
esis at joint probability > 99.73% CL. Only one multimask (multi- 
mask no. 53 of the HP 4 pixelization scheme) is used for this map, 
since only the one out of NrNm = 600 yields |P(Xq)l < Pth = 
1 — 0.9973. (See text and Table |IV]for more details), b) The spec- 
tral dependence of the kurtosis in the depicted regions in the three 
WMAP frequency bands along with the two and three sigma con- 
tours, from 1000 simulations, and the simulations mean are also plot- 
ted. 



In case of two other multi-masks, using which the INC data 
yield a detection at confidence levels of 99% and 98% in K 
(LB 64 8) and S (HP 8) respectively, the "3cr"-deviating re- 
gions turn out to be similarly located (like those in Fig.fTTb). 
With these regions masked out from the analysis a good con- 
sistency with the GRF simulations is reached. Note that the 
region 160 is rather small - mostly removed by the galactic 
sky cut (see Fig. \TTk and Table. |IV] for precise coordinates 
and size) however masking it has a comparable effect on the 
joint probability increase as masking out the two other and 
much larger regions. The consistency of the INC data with 
GRF simulations increases from 0.18% (without removal) to 
1.1% and 0.9% with regions 160 and 104+185 removed re- 
spectively. Individual removal of regions 185 and 104 only 
increases the level of consistency by < 0.5%. The simultane- 
ous removal of all three regions increases the consistency up 
to 12%. 

a. Dependence as a function of multi-mask To further 
test the robustness of this detection we have generated two 
other HP 4 pixelization scheme sets of Nm = 100 multi-mask 
each: one by simply choosing the three rotation angles with 
the prescriptions given earlier, and the other by focusing only 
on the region in the rotation angle parameter space within ±5° 
around the original orientation of the multi-mask leading to the 
3(7 detection. 

With the first set we obtained results yielding a joint prob- 
ability P(Xq) < 0-05 with 3 multi-masks, while in the sec- 
ond we find that 25% of multi-mask yield P(Xq) < 0.05, 
and 4% yield P(Xq) < 0.01, with the strongest detection 
P(Xq) = 0.0035, of which the Ucr map points to the same 
three regions as depicted in Fig. [TTb . We note that the re- 
ported regions (160 and 185) are located in directions towards 



which the strongest deviations in the individual region statis- 
tics (Fig.|7]) were found. 

b. Dependence as a function of frequency and resolution 
In Fig. [TTb we present the spectral dependence of kurtosis in 
the regions depicted in Fig. [TTb . While there is a non-trivial 
spectral dependence in regions 160 and 185, with opposite tilt 
- red and blue respectively - there is almost no spectral de- 
pendence in region 104. 

We also check the dependence on the S/N ratio in the se- 
lected regions. For this purpose we downgrade all datasets 
and simulations to resolution rig = 64, which effectively in- 
creases the S/N ratio per pixel by a factor of 8. We redo the 
multi-region analysis lowering the minimal region pixel num- 
ber threshold down to N^ix > 10 and find that the minimal 
probability per multi-mask (P (xg)) corresponds to a rejection 
threshold of 99.5% CL (Table |IV]). As seen in Table HYl the 
individual region response to the resolution change is strong 
only in case of region 160 while in the two other regions it is 
rather small. The result is robust under variations of region 
pixel number threshold and the number of simulations used to 
probe the underlying PDFs (7VsimPDF(x2) ^ {1000,5000}). 
Masking-out regions 104 and 185 reduces the anomaly to 
~ 96% CL and as expected in this case, removal of region 
160 has basically no impact on this value. 

c. Summary The non-trivial spectral dependence and 
close galactic-plane alignment in two of the three selected re- 
gions (160 and 185) suggests presence of some residual fore- 
ground anomalies. In case of the regions away from the galac- 
tic plane (104 and 185) since the local oddity is insensitive 
to the S/N ratio change it is also unlikely that an unknown 
instrumental noise fluke generates them. While we will re- 
turn to the overall statistical significance of these findings in 
Sect. IV CI we note that the positive KE seems to be incon- 
sistent with the extended foregrounds interpretation of these 
detections, according to the results from Sect.HV] Also the 3a 
detection of the multi-region analysis appears only in the INC 
data, but is clearly weaker in other single band maps. 

Given that this detection results from just one particular 
multi-mask and is selected from the lower-tail end of a whole 
distribution of equivalent measurements, it is inconclusive as 
regards indicating whether this detection is not just a fluke. 
Given that, we report in particular region 160, whose removal 
leads the overall significance to drop below 3a CL, as a ten- 
tative detection noting that more sophisticated local statistical 
analyses (see Sect. IVll) could be invoked to back these results 
up or refute them. 



2. Variance large scale distribution 

In Sect. IV Al we analyzed the large scale power distribu- 
tion as measured via multi-masks with regions of angular sizes 
ranging from 6° to 30° (Fig.[8]and Fig. [9]). 

In this section we focus on the joint multi-region analysis 
of the variance distribution. The corresponding results are il- 
lustrated in the second row of Fig. \TU\ The data remain in 
excellent consistency with the simulations. 

In Sect. lIVDl we found that the modulation amplitude (de- 
fined by Eq. [3]) of ^^^^^=1024 = 0.114 would be rejected at 
> 99.9% CL, and we argued that the modulation parameter 
^^max=40 = 0.114 would statistically be difficult to exclude 
at CL higher than 92%. Using our main set of the multi-masks 
(TableU) we fail to detect, in any of the data sets tested, any sta- 
tistically significant anomaly, such as the claimed hemispher- 
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Fig 12: Results of the "multi-region" analysis with the dipole-corrected V (V5) data. Joint probabilities i^(Xq) fEq. IB7I I for 600 multi- 
masks. The probabilities calculated with the spectral data maps from channels Q, V, W, INC are plotted with red, green, blue and black dots 
respectively. WMAP5 V5 data are plotted with green crosses (+) in HP 4 only. The V dataset original dipole has been replaced by the simulated 
dipole and removed in case of V5 data. Each dot (cross) corresponds to the joint probability of one multi-mask. From left to the right, the 
panels show the results with increasing resolution of pixelization schemes (see Table U) with 100 different multi-masks in each panel. Only the 
"mean" data is shown since all other results remain almost unchanged. For better visualization in range m G [10, 90] we plot only every lO'th 
sorted probability value. 
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Fig 13: Results of the multi-region analysis, using the set of 96 LB 
1 2 multi-masks, of the filtered in SH space up to imax= 40, low- 
resolution V data (left), and V5 data (right). The color in each 
pixel encode the multi-region, joint "probability of exceeding" de- 
rived with the LB 12 multi-mask rotated to the direction of the center 
of that pixel. Note that half of pixels in each map is redundant. 



ical power asymmetry (depicted in Fig. [8]), as measured from 
the large scale variance distributions in the multi-region anal- 
ysis. 

As an extension to that, we repeat the analysis performed 
in Sect. lIVDl for the low resolution, filtered in SH space up to 
^max= 40, WMAP data, using the same set of LB I 2 multi- 
masks, to test the variance distributions in the corresponding 
set of 96 differently oriented hemisphere pairs. We thereby ex- 
tend the test for the largest possible scales of 180°. We merge 
the Kp03 sky mask with the LB multi-masks for the analy- 
sis of the WMAP3 data, and KQ75 sky mask for analysis of 
the WMAP5 data. The result of the multi-region (here only 
two region) analysis is plotted in Fig. [13] for the V data (left) 
and V5 data (right). The minimal "probabilities of exceed- 
ing" found are: min (P(Xq)) ^ 3.3% (also marked in Fig.[6l 
bottom) towards b) = (247.5°, -30°)'^ for the V data and 
min(P(Xq)) ^6.5% towards (/, 6) = (281.5°, -19°) for the 
V5 data. 

Thes e two results agree well with the previously esti- 
mated teriksen et al.l l2007h intrinsic modulation parameter 
value A = 0.114 at scales ^max< 40, as they lie well within 
"one-sigma" region of the distribution of log-likelihoods ob- 
tained from 1000 simulations modulated with the modulation 



of A40 = 0.114. However we note that it will always be diffi- 
cult to reject such modulation at high confidence level as it is 
also realized (to this or a greater extent) on average in ~ 8% 
of GRF simulations (Sec. UVDl) . 

Note that the distribution of the probabilities of the joint 
multi-region analysis (Fig. [13]) has a very flat and extended 
maximum and, for example, the minimal joint-probability in 
the V data in the reported direction is only 0.2% smaller from 
the probability corresponding to the direction close to the 
galactic pole, which is roughly 50° away from the minimal 
probability direction. 

Analogous analysis, involving the LB I 2 multi-masks, 
but performed on the full resolution unfiltered WMAP V5 
data, results in larger minimal-probabilities: min (P(Xq)) ^ 
9.6% and the probability is minimized towards (/,6) = 
(225°, -78°). 



3. Residual dipole of the WMAP V channel. 

In Fig. [To] we see a deviation in the distributions of the 
mean in the V dataset (green dots) and in V5 dataset (green 
crosses) for most of the multi-masks. The fact that it is visible 
in most multi-masks suggests that the anomaly is not particu- 
larly sensitive to the multi-mask orientation, and that it comes 
from large angular scales. Indeed, as measured by the Ua- val- 
ues in individual regions of the multi-mask with the lowest 
joint probabilities, no region significantly deviates from simu- 
lations. 

However, we find that the V data are fully consistent with 
the GRF simulations if we remove the dipole from the data, 
which is roughly ~ 2 times larger than the one in our sim- 
ulations ^. Actually the dipole values in the datasets as 
measured by the multipoles / = 1 on the Kp03 cut sky 
power spectrum are: 47(/iK)^, 54(/iK)^, 45(/iK)^ in Q, 
V, W datasets respectively. The corresponding values in 
the WMAP5 data are: 64(/iK)2, 54(/iK)2 in V5 and W5 
maps respectively. The measurements of dipoles on ini- 
tially dipole-free maps, using a sky mask, introduce a bias 



^ This and the following result is accurate to within the tolerance of about ~ 
±7° resulting from the low-resolution search involving only 96 directions 
over a hemisphere in the HP 4 pixelization scheme. 



During the final stage of this w ork this sort of anomaly was independently 
reported in lEriksen et all i2008l) 
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Fig 14: Distribution of dipole orientations generated due to Kp03 cut 
sky from 1467 full sky simulations with vanishing initial dipole. The 
large spot in the center (at the right-hand side) of the plot, indicates 
the orientation of the WMAP3 (WMAP5) V band dipole. The color 
scale reflects the amplitudes of the leakage generated dipoles. 



due to power leakage from other coupled multipoles, lead- 
ing to non-zero dipole amplitudes. When sky-cut-generated 
dipoles are statistically accounted for, the result would yield: 
20ir(MK)^ 27t^r(38tr)(MK)^ 18ir(27ir)(MK)^ 
in the WMAP3 (WMAP5) data at 95% CL, and hence is con- 
sistent with vanishing intrinsic dipole (except for the V band 
channel). We note that the noise component generates dipoles 
with amplitudes of order Ci ~ 0.01 (/iK)^ which is about 
three orders of magnitude less than the leakage effect. How- 
ever the 95% CL effect is not sufficient to explain the strong 
anomalies detected in the regional tests. 

The anomaly is more visible in the difference of dipole 
amplitudes between different channels. The difference of 
9 ± 3(/iK)2 95^^ between channels V and W is ex- 
cluded using simulations at > 99.9% CL assuming that it is 
generated only by the power leakage from the cut sky. The 
difference in the WMAP V5 data is even larger: 11 ± 3(/iK)^. 

As for the amplitude of the V band dipole again, it becomes 
anomalous as one considers not only the magnitude, but also 
its orientation. The dipoles generated due to power leakage 
are strongly aligned within the galactic plane (Fig. [14]) as a 
result of the shape of the Kp03 sky mask. While the Q and 
W dipoles measured on the cut sky are close from each other 
and close to the galactic plane pointing at (/, h) = (13°, —8°) 
and = (7°, 5°) respectively, the dipole of the V band 
points at (/, h) = (350°, 30°) which is itself anomalous at > 
97% CL. WMAP5 data yield the dipole orientation (/, 6) = 
(203°, 28°) (Fig.m. 

We note that all dipoles with 6 > 30° have much smaller 
(roughly by an order of magnitude) amplitude than the one in 
the V band of the WMAP, which we believe is the reason for 
strong detections in the regional statistics carried out in the 
previous sections. When combining the alignment of the V 
band dipole with its magnitude, the hypothesis that it's gen- 
erated only via the power leakage can be excluded at a very 
high CL since out of 1300 simulations, and within the sub- 
set of 37 that have generated dipole aligned at |6| > 30° the 
maximal generated power is of only 13(/iK)^ which makes 
even the CL estimation unfeasible, since this when compared 
to the 54(/iK)^ of the V dataset (65(/iK)^ for V5), the simple 
test implies a rejection basically without doubt to a very 
reasonable limits. 

By reducing the dipole amplitude to the level consistent 
with simulations, or alternatively, by replacing it with our of 
our simulated dipoles, the data become consistent with our 
simulations at < 2a CL in the joint multi-region analysis (see 



Table V: Results of the dXX-multi-masks analysis for the signal dom- 
inated (co-added) and noise dominated (difference) maps. The 
columns content is as follows: (1) data set, (2) pixelization scheme, 
(3.. 5) "probability of rejecting" the consistency with GRF simula- 
tions (Sect. IB 31 ) for each MOD. In case of V and V5 datasets, the 
probabilities for the data with corrected dipole component are given 
in brackets. We round to integer percentiles for probabilities < 99%. 
The saturated values are marked with We abbreviate the results 
consistent at given CL as: "no significant detections (CL)". 
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Fig. [121) and at < la CL in all-multi-masks analysis at all 
resolutions (see Table |Vl in the next section). Note that the 
presence of this dipole in the V band is of no cosmological 
consequences since the dipole is marginalized over for any 
cosmological analysis, but may be important for other low-^ 
analyses. 



C. All-multi-masks analysis 

We now discuss the result of the all-multi-masks analysis 
described in Sect. IB 31 The corresponding results are pre- 
sented in Table El We see that the data are consistent with 
the simulations at ~ 68% CL. The previously mentioned (in 
Sect. IV B Tl) tentative NG detection in the INC data in kurtosis 
has indeed the largest probability of rejecting across the scales 
(61%) but it turns out to be statistically completely insignifi- 
cant. 

The large scale variance distribution is found to be perfectly 
consistent with the GRF simulations. 

We find a significant anomaly in the dipole component of 
the V band channel (see also Sec. IVB 3b detected via distri- 
bution of means in both WMAP3 (99.8% CL) and WMAP5 
(99.3% CL) data. 

All of the difference maps (see Sect. IV ED show a very 
strong departures from foregrounds free, simulated, difference 
maps, most prominently detected in means distributions. 



D. The "cold spot" context 

A NG anomalous kurtosis excess of a wavelet conv olution 
coefficients has been reported (e.g. ICruz et all (l2007h ) in the 
southern hemisphere, and was found to be associated with the 
locally cold spot (CS) in the CMB fluctuations around galac- 
tic coordinates (/, b) = (209°, -57°). In that work the wavelet 
convolution scales ranging from ~ 6.6° to ~ 13.2° in diam- 
eter were used, with anomaly being maximized at scales of 
~ 10° with a rejection on grounds of Gaussianity assumption 
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Fig 15: Kurtosis n-sigma map, thresholded at na,th = 3, combined 
from 100 multi-mask of the HP 8 pixelization scheme. The "cold 
spot" is marked with "CS". 



exceeding 99% CL. At the same time the authors note that the 
CS was not detected in the real space analyses. 

The range of scales mentioned correspond roughly to the 
scales tested by the HP 4 (- 14.6°) and HP 8 (- 7.3°) pix- 
elization schemes (see. Table U). 

In Fig. [15] we plot the scrambled n-sigma map of kurto- 
sis from the single-region analysis, obtained from 100 multi- 
masks of the HP 8 pixelization scheme, which most closely 
corresponds to the scales, at which the CS was detected. The 
deviation of — 3.6<j in the "the cold spot" direction, centered 
at (/, 6) = (209°, -53°) in Fig.[T5]is clearly found along with 
many other, locally significant, deviations. This particular CS 
however is not found at the same location in e.g. HP 4 pix- 
elization schemes or in skewness n-sigma maps, but rather it 
is shifted towards smaller galactic longitudes. 

As an extension to the tested scales, in this section we run 
a high resolution test using pixelization scheme HP 16 corre- 
sponding roughly to scales of ~ 3.7° and the INC data. We 
use 10 additionally generated multi-masks in this resolution, 
and we perform the single region, joint multi-region, and all- 
multi-masks statistics. We also performed the same analysis 
using the filtered up to ^max= 40 is SH space, low resolution 
(jis = 64) maps in which the spot is clearly visible. 

Although we find a locally negative KE and positive ex- 
cursions from expected distributions by > 3cr around the CS 
direction in variance, we also find similar excursions at several 
other directions. The CS itself is well localized in the scram- 
bled n-sigma maps of means with minimal value —2.9. 

However, none of these detections (see also Fig. [7]) hold un- 
der the scrutiny of the multi-region and all multi-masks analy- 
ses (Table El), that find these local anomalies to be statistically 
insignificant. 



E. Differential maps tests 

We discuss results of the QV, VW and QW difference maps 
tests as a simple cross-check with the CMB signal dominated 
maps tests, and a rough estimation of the residual foregrounds 
amplitude. We limit the number of multi-masks to 10 and 
use A^sim = 10^ simulations in single region analysis and 
A^simPDF(x2) = 10^ (A^sim(Cp) = 9 • 10^) simulations in joint 
multi-region analysis as before. 

As shown in Table |V] the residual foregrounds are very 
strongly (> 99.9% CL) detected in all difference maps, due 
to anomalies in means distributions. In particular the n-sigma 
and difference (A = m — (m)) maps of means distributions 



Table VI: Residual foregrounds amplitudes in the cross-band dif- 
ference maps. The columns contain: (1) pixelization scheme, (2) 
approximate angular size of the region as inferred from number 
of regions, (3) approximate corresponding multipole number ^ — 
ISO/Qreg, (4.. 6) - standard deviation of the difference maps outside 
Kp03 sky mask. 

Reg.sk. ^7reg [deg] i (j{qy) [/iK] (7(VW) [/iK] (7(QW) [/iK] 
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in QV, most prominently exhibit a dipole like structure ori- 
ented roughly at (/, 6) = ( 260°, 60°) which is close to the 
kinetic dipole direction (Fig. [T6h). The VW n^r map has a 
similar dipole structure, but with opposite orientation, which 
however is absent in the QW map. We find this to be a 
consequence of the previously detected (Sect. IVB 3b anoma- 
lous dipole component of the V band channel. Removing the 
dipole components from the difference maps, we have redone 
the three stages of the analysis, and although the dipole struc- 
ture ceased to dominate, we still find a very high (> 99.9% 
CL) rejection probabilities for in means as quoted in Table FVl 

In Table |Vll we present the amplitudes of the residual fore- 
grounds as measured by the variance of the A difference maps 
of means distribution for different scales as probed by our 
pixelization schemes. These remain in a good consistency 
with the limits given in lBennett et al.l (l2003h for residual fore- 
grounds contamination. 

In Fig. [T6tc-h) we show the rifj maps with distribution of 
the regions in the difference maps outstanding from simu- 
lations at significance larger than Sa (i.e. we use Ucj^th = 
3). Clearly, the close galactic plane regions are strongly de- 
tected. We note that the the KQ75 sky mask partially removes 
the most affected regions around (/,6) = (233°,— 10°), 
{l,h) = (259°, 18°) and the previously-mentioned (/, 6) = 
(199°, -23°). 

It is interesting to note that the largest scale negative 
(ricr < 0) anomalies seen in HP 2 (Fig. [T6b-d) away from 
the galactic plane, can also be induced by the foregrounds 
dominating along the galactic plane (with Ucr > 0) due to 
a very strong linearity of foregrounds induced quadrupoles 
with strong maximums aligned along the galactic plane 
and consequently strong minimums allocated close to the 
poles (Fig. [T6b). Such mechanism of foregrounds-generated 
linearity of the quadrupoles of the difference maps will not 
work in the foregrounds-free simulations, adding thereby to 
the observed large scale anomalies as probed via the largest 
regions. This effect is considerably smaller in the higher 
resolution pixelization schemes. 

In order to test the consistency of our noise simulations with 
the noise of the WMAP data, and the approximation the uncor- 
rected, white noise and to constrain limits of the systematical 
uncertainties, we also performed analysis using Q12 and VI 2 
difference maps in HP 2 pixelization scheme. The details of 
this analysis is given in Appendix O Here, we briefly report 
the result that the systematical effects measured, as before, by 
the standard deviation of the difference A maps remains at 
level < 1.7/iK at the scales corresponding the HP 2 pixeliza- 
tion scheme i.e. > 30°. 
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Fig 16: Results of the single-region difference maps analysis. 
In panel a) the residuals (Am(QV) = m — (m)) for the dif- 
ference QV data, uncorrected for the anomalous V band dipole 
is plotted. In panel b) the quadrupole of the VW difference 
map and in panels c) through h) for the three HP pixelization 
schemes, we plot the thresholded at 3cr, ria maps of means dis- 
tributions of the differential datasets QV (left) and VW (right) 
with the anomalous dipole component removed from the data. We 
make these maps (and maps of higher MODs) publicly available at 
|http : //cosmo . torun . pi/ ~blew/SKregstat /| 



VI. DISCUSSION 

A. Sensitivity, correlations and extensions 

Although we have shown that the statistics is rather help- 
less to robustly detect the NG signals considered (defined in 
Sec. IIVBI) via MODs higher than the variance, the statistics 
also proved to be a sensitive and precise tool for statistical 
isotropy measurements via variance and means. While we fail 
to detect the NG templates (inducing signals of rms ~ 100/iK 
within spots of ~ 10° via skewness and kurtosis) in the 
multi-region NG analysis, the single region analysis detects 
these as locally significant (Ua- ^ 3). Instead such template 
can be detected in the diW-multi-masks analysis at > 99.8% 
CL via variance. 



prominently seen in case of the reported dipole anomaly in the 
V band of the WMAP data (Fig. [T0|, where it is easily seen 
that depending on the choice of the pixelization scheme (e.g. 
such as those associated with the results at the right-hand side 
of the top plots in Fig. [T0| the obviously strong anomaly can 
be overlooked. 

It is important to mention the correlations between various 
multi-masks. This correlation occurs since, although the multi- 
masks sample the data differently, in the end, the same data 
are being sampled. The degree of redundancy is directly con- 
nected to the number of multi-masks used in the analysis and 
the magnitude of the correlation is related to the size of the 
regions used in a pixelization scheme. To quantify this we 
carried out a simple statistical test only whose goal was to es- 
tablish whether our test is statistically more sensitive to the 
change of the tested data itself or to the change of the multi- 
masks for various resolutions parameter r and MOD. 

If the results from different multi-masks were strongly cor- 
related, then for a constant number of DOF one would ex- 
pect the variance of the results measured over these multi- 
masks (e.g. values) to be much smaller than the variance 
computed while fixing a multi-mask but varying dataset. On 
the other hand, if the variance of the results when changing 
dataset was much smaller than the variance when changing 
multi-mask, then the test would not be very sensitive to par- 
ticular features in the data, or even unstable. To test this we 
calculated the R statistic defined as follows: 

(cTsim (x^ {r, MOD) /DOFeff ) ) 

i?(r, MOD) = ; ' (4) 

^ ^ (a^(x2(r,M0D)/D0Feff))^^^ 

where () . denotes an average over simulations and () 

\/ Sim ^ \/m 

denotes an average over multi-mask for a given simulation, 
and where DOFeff = DOF{r, m) is the effective number 
of degrees of freedom. Measuring this R statistic using our 
simulations, we find that the test is approximately equally 
sensitive at all resolutions and for all MODs, and that eventu- 
ally it is little more sensitive to the change of the data under 
test than to the change of the multi-mask giving R values 
around 1.2. 

The approach with arbitrary shape of the regions {multi- 
mask) and their orientation is quite flexible, and different 
shapes can possibly be used for different applications inde- 
pendently of the enforced sky cuts. This allows for a thorough 
test of the multivariate nature of this Gaussian field. One could 
also consider other statistics than the first MODs, as e.g. re- 
gional Minkowski functionals. Another possible extension is 
to apply a specific pre-filtering of the data in the SH space in 
order to expose for the test features dominating at particular 
scale. Such slicing of the data into subsets of maps according 
to some chosen ranges of multipoles could in principle sig- 
nificantly improve the sensitivity. The multi-region full- sky 
analysis though is restricted generally to the resolutions up to 
which the full co variance matrix analysis is feasible. 



In lEriksen et al.l (|20Q4|) a regional statistical analysis was 
performed using a set of circular regions uniformly distributed 
across the sky. Our analysis is similar in spirit but uses differ- 
ent statistics and a richer sets of regions, varying both in size 
and shape. This approach has been validated by the fact that 
we have shown that the resulting statistical signal can be a sen- 
sitive function of the particular choice of regions. This is most 



VII. CONCLUSIONS 

We introduce and perform a regional, real space test of sta- 
tistical isotropy and Gaussianity of the WMAP CMB data, us- 
ing a one-point statistics. We use a set of regions of vary- 
ing size and shape (which we call multi-mask) allowing for 
an original sampling of the data. For each of the regions we 
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analyze independently or simultaneously, the first four mo- 
ments of distribution of pixel temperatures (i.e. mean, vari- 
ance, skewness and kurtosis). 

We assess the significance of our measurements in three 
different steps. First we look at each region independently. 
Then we consider a joint multi-region analysis to take into ac- 
count the spatial correlations between different regions. Fi- 
nally we consider an ''all-multi-masks'' analysis to assess the 
overall significance of the results obtained from different sky 
pixelizations. 

We show that the results of such multi-region analyses 
strongly depend on the way in which the sky is partitioned 
into regions for the subsequent statistics and that our approach 
offers a richer sampling of the underlying data content. Conse- 
quently the diW-multi-masks analysis provides a more robust re- 
sults, avoiding possible biases resulting from an analysis con- 
strained only to a single choice of pixelization scheme. 

We find the three-year WMAP maps well consistent with 
the realistic, isotropic, Gaussian Monte-Carlo simulations as 
probed by regions of angular sizes ranging from 6° to 30° at 
68% confidence level. 

We report a strong, anomalous (99.8% CL) dipole "excess" 
in the V band of the three-year WMAP data and also in the V 
band of the WMAP five-year data (99.3% CL) (Sect.lVBSj). 

We test the sensitivity of the method to detect particular 
CMB modulation signals defined via the scale dependent mod- 
ulation amplitude parameter {Ai^^^) for the case of modula- 
tion extending up to maximal multipole number of £max = 40 
and £max = 1024. We are able to reject the modulation of 
amplitude of A1024 = 0.114 at > 99.9% CL and find that 
A40 = 0.114 can be statistically excluded only at ~ 92% CL 
(Sect. HVDl IVB2l) . Given the WMAP V band data, we find 
that the large-scale hemispherical asymmetry is not highly sta- 
tistically significant in the three-year data (~ 97%) nor in the 
five-year data (~93.5%)at scales ^ < 40. Including a higher- 
l multipoles only decreases the significance of hemispherical 
variance distribution asymmetry. 



We also test the sensitivity to detect a broad range of small 
(10° in radius) locally introduced NG signals, inducing non- 
vanishing kurtosis (and in general skewness) of rms amplitude 
~ 100/iK and find that the method is able to detect these as 
locally significant, but the overall impact in the joint multi- 
region analysis is unnoticed by mean, skewness and kurtosis, 
but is strongly detected (~ 99.8% CL) by variance distribu- 
tions. We conclude that the NG foreground-like signals will 
be easier to detect using local variance measurements rather 
than higher moments-of-distribution. 

We also analyze our results in context of the significance of 
the "cold spot" (CS), reported as highly anomalous at scales 
corresponding to ~ 10° in diameter. While we notice the cold 
spot region as having locally anomalous, negative kurtosis- 
excess and locally increased variance (eg. Figs. [71 [15]), we do 
not find these deviations to be globally statistically significant. 

We easily detect the residual foregrounds in cross-band dif- 
ference maps at average rms level < 7/iK (at scales ^6°) 
and limit the systematical uncertainties to < 1.7/iK (at scales 
> 30°) as a result of the analysis of same-frequency differ- 
ence maps. These levels are consistent with the previously 
estimated limits. 
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Appendix A: WMAP SIMULATIONS 
1. Signal and noise pseudo tests. 
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Fig 17: Consistency check between the simulations and WMAP 
three-year observations. The pseudo- power spectra of the WMAP 
(light blue/dot-dashed curve) and its simulation (the underlying 
blue/dashed line) in channel Ql are well consistent both in Kp03 sky 
cut regime and in the noise dominate d regime. Reconstructed power 
spectrum of the lHinshaw et al.l (l2007l) is plotted as big red crosses. 
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In order to assess the confidence levels we have per- 
formed Monte-Carlo simulations of the signal {^max < 1024 
and Co,i = 0) and inhomogeneous but uncorrected Gaus- 
sian noise maps for all DAs, according to the best fit to 
the AC DM model power sp ectrum, extracted from obser- 
vations (lHinshawetal.ll2007l) . The simulations include the 
WMAP window functions for each DA. The A/'sim(= 10"^) 
full sky simulations were generated at the HEALPIX resolu- 
tion Us = 512, for each DA and preprocessed in the same way 
as the data. 

In Fig.[T7]an example of simulated, and recovered pseudo- 
Ci is compared with the pseudo- Q of the WMAP data from 
channel Ql as a simple consistency check. Similar simulations 
were performed for the remaining DAs. 

As a simple consistency test, we present in Fig. [TSl statistics 
of the full sky WMAP third year data as compared with the 
simulations. The data are well consistent with the simulations. 



a. Tests of modulated- simulations ' power spectrum 

As mentioned in Section UVDI the modulation (Eq.O alters 
the underlying power spectrum of the modulated simulations 
and could possibly mislead us in the interpretation of the high 
rejection confidence level thresholds, reported in Table [III for 
the modulated simulations, since the additional power could 
be merely a result of the inconsistency on grounds of the 



Fig 18: The distributions of means, variances, skewness and kur- 
tosis of A^sim = 10^, full sky, INC simulated data realizations, cal- 
culated outside the Kp03 sky mask. Skewness and kurtosis values 
of the distributions are also given. Vertical bars indicate the values 
derived from the WMAP three-year data. The quantile probabili- 
ties of the mean, variance, S, and K values of the WMAP data are 
{0.57, 0.015, 0.34, 0.35} respectively - well consistent with Gaus- 
sian, random simulations. Low probability of the total variance re- 
sults primarily from the low quadrupole (octupole) of the WMAP 
data as compared with the best fit AC DM model. Note that the 
distribution of the means of the simulations represents only a numer- 
ical noise since during preprocessing all maps were shifted to zero 
the mean outside the Kp03 sky mask (< T >^ 0) and therefore it 
does not carry any important cosmological information. The spectral 
analysis yields a similar results. 



intrinsic power spectrum mismatch, rather than the regional 
variance distribution analysis, and violation of the statistical 
isotropy. 

In this section we quantify this effect. We gen- 
erate a set of 10 WMAP V3 simulations and mod- 
ulate their CMB component using modulations 
Ai024(n = (225°, -27°)) = {0.114,0.2}. We next 
extract their Kp03 cut sky pseudo-power spectra up to 
^max = 700 i.e. up to the scales where the noise compo- 
nent already strongly dominates over any possible CMB 
modulation signals. We derive the contribution to the total 
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variance of the map per multipole according to: cr^ = ^^^Ci 
(£ e {2, . . . , 700}). We measure the degree of the consistency 
of the modulated simulations' power spectra with the non- 
modulated simulations' power spectra, using the unbiased 
estimator of the full covariance matrix E/^/ = Cov(cr^, a^/) 
derived from 3000, analogously prepared a£ vectors extracted 
from the GRF WMAP V3 simulations, and using the corre- 
sponding Monte-Carlo probed values distribution from 
223 additional simulations. 

We find that the average rejection thresholds for the mod- 
ulations Aio24 = 0.114 and A1024 = 0.2 are 54% and 59% 
respectively with the standard deviation ~ 30% in both cases. 
We therefore conclude that our results given in Table HIl can- 
not result from simply alteration of the power spectra after the 
modulation field has been introduced. Similar results are ob- 
tained if the off-diagonal terms of the covariance matrix are 
neglected, which indicates that the low detection threshold 
does not result from the limited number of simulations and 
a possible non-convergence of the covariance matrix. 



Appendix B: ASSESSING STATISTICAL SIGNIFICANCE 

In this appendix we give details on out statistical ap- 
proaches. For fast reference, we summarize the principal sym- 
bols used in Table: IVIII 



1. Individual region analysis 

In the case of individual regions statistics for every ith 
region (i G {1^ ..^ Nreg{r^m)}) of every multi-mask m G 
{l..A^^} and for every MOD (X G {m, a, S, K}) and every 
dataset {d G {Q, V, W, INC}) we define a parameter vector 
p = {r, m, d} and independently calculate the tail probabil- 
ity P(Xp ^) of occurrence of Xp ^ value of the data in the 
^sim = 10^ simulations, probing the corresponding probabil- 
ity distribution functions (PDF). The quantities m, a, S, K cor- 
respond to the first four moments of distribution respectively. 



Hence we define P(Xp^i) as: 

p(Xp,o ^ p,{\x;j - Q2p,| > \x^T - Q2p.iI) (Bl) 

where (52p,i is the second quartile of the corresponding PDF, 
and Pq is the quantile tail probability. 

In principle, considering N simulations allows us to probe a 
region of the underlying PDF corresponding to Gaussian num- 
ber of "sigmas" 

± n^^ = V2erf-^(1 - 2/N) = |cdfG"'(l/N)| (B2) 

where cdfc"^ is the inverse Gaussian cumulative distribution 
function (CDF). For N = Nsim = 10^, as it is the case for 
individual region statistics, ^3.72 corresponding to a 

probability of exceeding of 0.02%. 

To derive the probability from Eq. lBll we use linearly inter- 
polated quantile probability within the MC probed PDF range: 



Qiin{Xp,i) for 

(min(X^^),max(Xj 
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where Qun is the linearly interpolated quantile probability, 
and min(Xp^) and max(Xp^) are the minimal and maxi- 
mal Xp i values obtained from a sample of N^im simulations 
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Fig 19: Example of extrapolation (solid lines) used to compute the 
distribution of 10^ MODs from one of the regions in one of the multi- 
masks. Mean and variance MODs values were multiplied by a factor 
of 10^ for the sake of clarity. Note that the Gaussian extrapolation 
(used only outside MC probed PDF region) also interpolates the data 
quite well in case of the means (red crosses). 



that probe the underlying Xp^^ PDF. Outside the probed range 

(Xp,i G (-oo,min(X^^^)] U [max (X^^^), 00)) we extrapo- 
late using a Gaussian distribution of the form 



P X, 



= (l-erf(i^ 
n^C(l 



V2^ 



(1 



|min(X--)-Q2p,J^ 
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> max(X;^7) 
(B4) 



is connected to the 
matching condition 



Note that the extrapolation form 
MC probed PDF range by the 

P(Xp,i(n^^)) = P(min(X^;j^)) = P(max(X^^)) = 2/N 
- i.e. the probability of finding a Xp,, value anywhere in 
range (—00, min(Xp^^)] U [max(X^), 00). An example of 
this extrapolation is shown in Fig. Q9\ This extrapolation is 
obviously not validated for MODs higher than the mean, but 
we use it as a guide for very low probability events. Note that 
all our results with a lower significance (roughly outside 3a 
CL) are obtained modulo this approximation. 



2. Multi-region joint analysis 

In the multi-region analysis, we account for all correlations 
between regions of a given multi-mask using an unbiased esti- 
mator of the full covariance matrix. Using the same param- 
eter vector p = {r, m,(i}, we define a one column vector 
for each MOD (Xp G {mp, cTp, Sp, Kp}) of size Nreg{r^ m) 

such that Xp = (X^,^,(^,^^i, X^,^,c/,i=Ar^eg(r,m)) con- 
tain Xp ^ values for each region of a given multi-mask m of 
given pixelization scheme r and dataset d e {Q, V, W, INC}. 
Introducing a parameter vector q = {X, r, m, (i} we define a 



corresponding Xq value as: 



Xq 



(x: 



data 
P 



Sim 
P 



))^c-i(x^-*- 



(X;™)) (B5) 



where the () is the average X from 



sim(Cq 



simulations and 



Cq -"^ is an unbiased estimator of the inverse covariance ma- 
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Fig 20: Convergence of Xq values as a function of the number of sim- 
ulation used for the covariance matrix calculation for a given MOD. 
Each color corresponds to a type of pixelization schemes. The num- 
ber of regions increases from bottom (48, for HP 2) to top (1024 
for LB 64 16) (see. Table U). Horizontal lines indicate the effective 
numbers of degrees of freedom for that multi-mask (i.e. number of 
unmasked regions due to Kp03 sky mask). These corresponds to 
the the first MOD only (mean) but other MODs exhibit similar de- 
pendence. 



trix (iHartlap et al.l2QQ7h calculated from A/'sim(Cq) simulations 
and is given by: 



-^sim(Cq) - Nreg{r,m) - 2 - ^ 



sim(Cq 



and Cp ^ is the inverse covariance matrix. 

In Fig:[20lwe show the convergence Xq with the number of 
simulations used to calculate the covariance matrix, A/'sim(Cq) ^ 
for all six types of pixelization schemes. The biasing of the Xq 
values with regards to the ideal Xq(^sim(Cq) = oo) value is 
caused by the limited number of available simulations. As ex- 
pected, for a given number of simulations, this bias decreases 
with the effective number of DOF and as such with the multi- 
mask resolution. We account for this biasing by using sim- 
ulations to probe the underlying PDF of the Xq values, in- 
stead of using theoretical x^ distributions. In fact, it would 
not be valid to use theoretical x^ PDFs, since the distribu- 
tions of individual MODs (except for the mean) are not Gaus- 
sian. As such, we will probe the underlying distributions using 



sim(Cq 



simulations, for each MOD 



simPDF(x^) — ^'^sim 

and for each multi-mask. 

In the case of the joint multi-region statistics we use N = 
A^simPDF(x2) = 10^ simulations corresponding to ^ 
3.09, which gives a probability of exceeding and 0.2%. The 
remaining A/'sim(Cq) = 9 • 10^ simulations are used for covari- 
ance matrix calculation. 

Given a (Xq)^^*^ value we define a corresponding joint 
probability as: 

p(x^)^i-p,((x^f->(x^)^^*^) 



(B7) 



We use the same formulas for inter/extrapolation as described 
in Sect. IB Il by replacing Xp i with Xq- 

We note that, in fact, in case of the multi-region analysis 
it doesn't make any difference whether the derived x^ values 
are statistically debiased or not, since exactly the same bias- 
ing affects the simulated x^ values used to probe the under- 
lying PDF. We also find that the convergence of such derived 



Table VII: List of principal acronyms used in this section and in 
the main body of the paper, briefly summarized for quick reference. 
^ - indicates "unless specified otherwise". 



Symbol 



explanation 



MOD 

data 
sim 



Nregir, TJl) 



P 

q 



MC 



n 
Xp 



sim(Cq) 



^simPDF(x2) 



(B6) P{xl) 



moment of distribution 

upper index to indicate a measurement on data 

upper index to indicate a measurement on simulations 

total number of multi-masks in a given 

pixelization scheme (Nm = 100*) 

multi-mask index number m G {1, . . . , Nm} 

pixelization scheme resolution parameter 

r e {48, 192, 256, 512, 768, 1024} (Tab.lB 

MOD:X G {m,cr, S,K}, 

mean, variance, skewness, kurtosis respectively 

dataset: d e {Q, V,W,INC} 

number of regions in m'th multi-mask 

of the pixelization scheme r 

region index within a multi-mask 

parameter vector: p = {r, m, d} 

parameter vector: q = {X, r, m, d} 

value of a MOD for i'th region of m'th multi-mask 

of r'th pixelization scheme in a dataset d 

total number of simulations; 

number of simulations used in the single-region analysis 
probability corresponding to Xp i derived from A^sim 
simulations (Eg. lBlt 
defined in Ea. lB2l 

vector of MOD values for for a given value of parameter p 
number of GRF simulations used to derive the covariance 
matrix in multi-region analysis (A^sim(Cq) = 9 000*) 
number of GRF simulations used to probe the distribution of 
Xq values in multi-region analysis (A^simPDF(x2 ) ~ ^ 000'^) 

value for the corresponding Xp^*^ (Eq. lB5b 
probability corresponding to Xq (Eq.|B7} 



probabilities is actually much better than one could infer when 
looking at the worst case of LB 64 16 presented in Fig.[20l We 
estimate that the derived probabilities converge to their true 
values within ~ 10% of that value in this worst case of LB 64 
16 pixelization scheme, while in case of the LLP 2 the conver- 
gence of the derived probabilities is < 2%. 
The statistical debasing used for the values matters how- 
ever in our third-stage analysis, i.e. in case of the dXl-multi- 
masks analysis. 



3. All Multi-masks analysis 

There is no a unique way to generalize from the results of 
the multi-region analysis. Although the most straightforward 
way would be to compute the inverse covariance matrix be- 
tween all the MODs and for all regions of all multi-masks, this 
turns out to be computationally not feasible. 

Note the fact that the x\ values are partially correlated with 
each other, since they result from a different sampling of the 
same dataset. However the degree of correlation strongly de- 
pends on the particular multi-mask properties and resolutions. 
In particular the mitr-multi-mask correlations are largest in 
the lowest resolution multi-masks. The smaller correlations 
between various multi-masks, the more additional information 
the multi-mask analysis explores about the dataset. Large 
intQY-multi-mask correlations indicate that not much new 
information is gained making it unnecessary to process large 
number of multi-masks. 

In the following, in order the integrate all the information 
probed by different multi-masks we calculate the cumulative 
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probability of rejecting the GRF hypothesis using the median 
distribution ^(Xq^^) of all processed Xq^c distributions and 
the median Xq value of the data. Therefore we calculate the 
distribution and Xq value as: 

^(4mc) = (¥'(X^Mc/DOFeff))„» 

_ (B8) 

Xl = {x^/DOFeff)„. 

where ()m is the average over all multi-masks, DOFeff = 
DOF(r, m) is the effective number of degrees of freedom 
which is a function of resolution r, of the pixelization scheme 
and of a particular multi-mask due to interplay with the sky 
cut. 

We define the joint cumulative probability of rejecting the 
GRF hypothesis of the WMAP CMB data via inconsistency 
with our simulations as a function of q analogically as in 
Eq. IB71 and we use the same extrapolation and interpolation 
formula as in case of multi-region analysis. 

Appendix C: NOISE SIMULATIONS TESTS 

Difference maps obtained from observations in nearly the 
same frequency, and with nearly the same beams profile pro- 
vide a good opportunity to measure the statistical properties of 
the instrumental noise. 

We have performed a reduced x^ tests, directly in pixel do- 
main, of the difference maps obtained from different channels 
of the WMAP data (Q12, V12, QV and QlVl) at the Healpix 
resolution rig = 512, and compared with results of the same 
tests performed at a low Healpix resolution Ug = 4. 

Since the covariance matrix of the noise realizations is well 
diagonal a single variate Gaussian statistics was assumed, and 
reduced x^ distributions used. The Q12 and V12 yielded 
a well consistency with the simulations at both resolutions. 
The QV and QlVl difference data however, turns out to be 
more troublesome. Whereas there is a good consistency at 



high resolution, the low resolution reduced x^ t^sts show 
significant discrepancy yielding a "probability of rejecting" 
P = 0.999963 in case of QV map and P = 0.998 in case 
of QlVl map. This result is also discussed in Sect. IV El in 
light of the anomalous dipole component of the V band map. 

We also performed a single-region, joint multi-region, and 
all multi-masks analyses on the Q12 difference map, using 
a subset of 10 selected multi-masks of the HP 2 pixelization 
scheme. Since the low resolution analysis yields a quick con- 
vergence (Fig. [20lin Sect. IB 21) a small number of 500 simu- 
lations were generated and half of them used for covariance 
matrix calculation, and the other half was used for probing 
distributions of the x^ values. 

We found strong anomalies in the distribution of means (of 
which joint probability is well extrapolated using Eq.|B4] out- 
side the MC probed range (Fig. [19])). The variance of the 
scrambled A maps show that the rms amplitude of the dif- 
ferences is limited to the 1.7/iK at these scales which is con- 
sistent with the limits to the residual sy stematical uncertaint ies 
in Ql and Q2 channels of the WMAP dHinshaw et al.ll2003h at 
these scales. The constraint includes not only the systematical 
uncertainties but also possible differences due to uncorrected 
white noise used in our simulations, which in principle in the 
regional statistics do not average out in the same way as the 
pre- whitened 1/f pink noise of the WMAP data. 

A difference A map of the variances can also serve as 
a rough estimate of the level of local systematical effects. 
Anomalies in the scrambled map are indeed found, with 
strongest deviations concentrated in parts of regions adjacent 
to the Galactic Center, with extreme values < 3/iK. How- 
ever a close orientation of the regions to the Galactic Center is 
more likely a hint on the residual foregrounds contamination, 
due to slight differences in the effective frequency of the Ql 
and Q2 differential assemblies, as well as in the beam profiles, 
rather than a manifestation of a systematical effects. Due to 
this leakege the limits to the aforementioned systematical ef- 
fects at the level of 1.7/iK should be considered as an upper 
limits. 



